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The significant orbital eccentricities of most giant extrasolar planets may 
have their origin in the gravitational dynamics of initially unstable multiple 
planet systems. In this work, we explore the dynamics of two close planets 



merical scattering experiments. We derive a criterion for two equal mass 



on inclined orbits through both analytical techniques and extensive nu- 
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planets on circular inclined orbits to achieve Hill stability, and conclude 
that significant radial migration and eccentricity pumping of both planets 
occurs predominantly by 2:1 and 5:3 mean motion resonant interactions. 
Using Laplace-Lagrange secular theory, we obtain analytical secular so- 



lutions for the orbital inclinations and longitudes of ascending nodes, and 
use those solutions to distinguish between the secular and resonant dynam- 
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ics which arise in numerical simulations. We also illustrate how encounter 
maps, typically used to trace the motion of massless particles, may be mod- 
ified to reproduce the gross instability seen by the numerical integrations. 
Such a correlation suggests promising future use of such maps to model 
the dynamics of more coplanar massive planet systems. 
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1 Introduction 



Over the last decade, our catalog of planets has increased twelve- fold. Over 120 extrasolar 
planets around main sequence stars in over 100 systems have been discovered, 1 of which at least 10 
such systems contain multiple planets. Only three of the planets around solar-type stars have been 
discovered by transit (Konacki et al. 2003; Konacki et al. 2004; Bouchy et al. 2004), and for only 
one other has a transit event been detected 1 (Charbonneau et al. 2000; Henry et al. 2000). All 
other of these "Extrasolar Giant Planets" (EGPs) have been detected by the Doppler radial velocity 
technique, which provides information about an exoplanet's semimajor axis and eccentricity, but 
not its inclination with respect to either other planets or the stellar rotation axis. The dearth of 
data for this one parameter, in addition with the near-coplanarity of the major planets in our Solar 
System, has directed most studies of multiple-exoplanet systems to neglect any nonzero mutual 
relative inclination of the planets when discussing the system's past, present, or future dynamical 
evolution and resulting stability (an exception is the dynamical study of the 47 Ursae Majoris 
system by Laughlin et al. 2002). 

The prospects for two or more planets emerging from a disk on inclined orbits have not yet been 
fully explored. Even in the simplest scenario - in which the planetary disk itself is flat - excitation 
of planetary inclination (and eccentricity) can be driven by mean-motion resonances with the gas 
disk (e.g. Goldreich and Tremaine 1980), although those are opposed by dissipative secular effects 
(Lubow and Ogilvie 2001). Recently, Thomas and Lissauer (2003) have used three-dimensional 
numerical simulations to argue that non-coplanar planetary systems could be common, provided 
that damping of eccentricity by the disk is relatively weak. The eccentricity evolution of planets 
embedded within disks is, itself, rather uncertain (Papaloizou et al. 2001; Ogilvie and Lubow 
2003; Goldreich and Sari 2003), and probably depends upon unknown properties of the dynamics 
of eccentric disks (Kato 1983; Lyubarskij et al. 1994; Ogilvie 2001). In less standard models - for 
example those in which the disk itself is warped (Terquem and Bertout 1996), or the planetary 
system is initially crowded (Papaloizou and Terquem 2001; Adams and Laughlin 2003; Nagasawa 
et al. 2003) - high relative inclinations between a small number (often two) of surviving planets are 
even more probable. Finally, Yu and Tremaine (2001) conclude that when an inner planetesimal 
disk is driven into a star, the runaway multiple planet system may not be coplanar. 

The wide use of the coplanarity assumption has also affected the types of mean motion reso- 
nances analyzed for EGPs. Resonances may be classified as purely eccentric, purely vertical, or a 
combination of the two. All three types play a crucial role in shaping the dynamical evolution of 
ring particles, asteroids, and Kuiper Belt Objects, but for larger bodies like satellites and planets 
eccentric resonances dominate; the 4:2 Mimas- Tethys resonance is the only known purely incli- 
nation resonance among satellites or planets in our Solar System (Peale 1999; Champenois and 
Vienne 1999a). However, Henrard et al. (1995) and Wisdom (1987) have shown that one must in- 
clude inclination effects when attempting to characterize resonances in planetary systems. Further, 
most resonance studies have been conducted under the guise of the restricted three-body problem 
(Peale 1986; Murray and Dermott 1999), or a variation thereof, due to its analytical tractability 
and practical application to dynamical situations in our Solar System. The resonance between two 
EGPs cannot be modeled as in the restricted three-body problem because of the unknown relative 

1 From the on-line Extrasolar Planets Encyclopedia, at |http: //cfa- www. harvard. echi / planets/catalog. html (Schneider 2004) 
as of June 11th, 2004. 
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inclinations between the EGPs and because neither mass can be treated as negligible. 

The coplanarity or small-inclination assumption has affected other types of dynamical analysis 
as well. Gladman (1993) produced pioneering work on establishing the criteria for stability of two- 
planet systems, but did not include nonzero planetary inclinations in his analysis. Traditionally 
the expansion of the disturbing function has been obtained by expanding about inclination values 
of 0° (Kaula 1962; Murray and Dermott 1999). Although some forms of the disturbing function 
expanded about high eccentricity values exist and have been applied to nonplanar problems (Roig 
et al. 1998; Beague and Michtchenko 2003), to our knowledge high inclination expansions have not 
yet been derived. Encounter maps have been used as a technique to trace the orbits of planetesimal 
swarms (Namouni et al. 1996) and outer solar system objects (Duncan et al. 1989), but have not 
yet been used to model highly inclined systems nor systems with two massive planets. 

We address several of the above issues in this work by presenting analytical and numerical results 
for two massive planets on arbitrarily inclined orbits. Section 2 provides a description of the Hill 
stability limit, and illustrates an extension to the inclined case. This extension forms the basis for 
the ~ 10 4 gravitational scattering experiments we performed, the results of which are presented in 
Section 3. In Section 4, we derive compact expressions for the secular evolution of inclined planets, 
and compare those to the simulation results to aid in analyzing the latter. In Section 5 we illustrate 
how Namouni et al. 's (1996) encounter map may be extended to two massive planets on inclined 
orbits. Section 6 discusses the implications of this work, and Section 7 briefly summarizes the 
results. 

2 The Hill Stability Limit 

2.1 Background 

Due to the computational expense and large phase space involved, a numerical study of even two- 
planet systems requires a carefully chosen set of initial conditions. We determined what conditions 
are most suitable for study by first considering the types of stability a system may have. As 
summarized by Gladman (1993), a system that is Lagrange stable is one where the planetary orbits 
are bound for all time. A system that is Hill stable is one where the orbits of planets cannot cross. 
Only Hill stability, however, is mathematically tractable. 

Marchal and Bozis (1982) applied the mathematical results of Hill stability to two-planet systems 
with a central mass much greater than the planetary masses. Gladman (1993) analyzed the problem 
further, providing explicit expressions in terms of orbital elements which approximate the minimum 
initial semimajor axis separation of two planets necessary to ensure Hill stability. We expand on 
Gladman's (1993) results by deriving a general expression for equal-mass planets and an arbitrary 
initial inclination. 

2.2 Setup 

The Hill stability limit may be considered in terms of the fixed points of the general three-body 
problem. These fixed points are the locations in space that separate qualitatively different dynamical 
behaviors of the planets. Typically, the fixed points are analyzed for the circular restricted three- 
body problem, in which the secondary mass is much greater than the tertiary. However, not as 
well studied are the fixed points of Hill's problem (Hill 1878; Henon and Petit 1986), in which the 
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secondary and tertiary masses are comparable. Hill's problem has wide application to solar system 
dynamics, and especially to two-planet systems. 

Consider the case of a two-planet system, with planets of mass mi and m 2 orbiting a massive 
central object m 3 (so that m 3 ^> m 2 and m 3 ^> mi). Then for mi > m 2 , Marchal and Bozis (1982) 
showed that the Hill stability limit is given by the solution to: 



4 mim 2 mim2 (llmi + 7m 2 

+ m 2 J 

2 (mi + m 2 + m 3 



1 + 3 '^7^ — 1 2^ + --- 

(mi+m 2 ) 3 3m 3 (mi+m 2 ) 

(2.1J 



rC /l, 



G 2 (mim 2 + miffi3 + m 2 m 3 ) 

where c is the angular momentum and /i the total energy of the system. Note that the left-hand-side 
(LHS) of this equation is an approximate expansion, which can be continued if greater accuracy is 
required. Since c and h are functions of the orbital elements, this equation is a (quartic) equation 
for the separation needed for Hill stability. 

2.3 The Circular Coplanar Limit 
2.3.1 Units 

Gladman (1993) provided an analytical solution for the critical separation in terms of the masses 
alone. We adopt his useful notation in this work by setting the following parameters and defining 
the following auxiliary variables: Let m 3 be the "central mass" so that m 3 3> mi and m 3 ^> m 2 . 
Gladman (1993) defined: 

Em\ m 2 
rrii = 1, Hi = — , /i 2 = — , and a = /ii + /i 2 . (2.2) 
m 3 m 3 



eii = 1, a 2 = 1 + A, and 5 = \/l + A. (2.3) 



7i = yjl-el, and G = 1. (2.4) 

The quantity A represents the separation of both planets and will be used throughout the rest 
of this work. We work in units where the sum of the masses and the gravitational constant are 
separately equal to unity. Since we ignore physical collisions between planets, the problem remains 
scale-free, and ai may be scaled straightforwardly to arbitrary values. For generality, henceforth we 
quote semimajor axis in these scaled units (i.e. ai = 1 initially), and times in units of the orbital 
period at a = 1. Because m 3 is considered to be the central mass, and through the choice of units 
and approximations, one may eliminate mi, m 2 , and m 3 in Eq. (j2.1|) in favor of /ii, /i 2 , and a. 

2.3.2 Statement of critical separation 

With Gladman's (1993) variables, the LHS of Eq. (j2.1|) can now be written as: 

1 + 3 |/^_ ^ 2 (llMi + 7 M2 ) + 
ai 3a 2 
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In route to rewriting the RHS of Eq. (|2.1j) by using the definitions in Eqs. (j2.2j) - (J2.4j) . one must 
first re-express the total energy and angular momentum of the system as such: 



h " -y - 2<tta) • (2 - 6) 



Gm 1 ym 3 a 1 (1 - ef) + Gm 2 \/ m 3 a 2 (1 — e|) + • • 
A*i7i + ^272(5- 



(2.7) 



Now the RHS of Eq. ()2.1|) can be written as: 

a ' 3 fa + ff) O^i + A^) 2 . (2.8) 

Equating ()2.5|) and (|2.8|) for circular orbits yields the critical separation of two planets in terms 
of the planet masses only (Gladman 1993): 



A crit = 2-36(^1 + ^2)3 +2 



o±/ \| ll/il + 7/!2 



+ ■ ■ • (2.9) 



3 e (/ii + /i 2 ) 

If two planets are initially separated by a distances greater than A cr i t , then their orbits cannot 
cross for all time, and the planets are said to be Hill stable. In order to evaluate Eq. ()2.9|) for 
typical giant planet mass ratios of \i — ~ 10~ 3 , the equation may be rewritten in a more practical 
way, as: 

^^Mi^ + i^)'^ 1 ) (2 ' 10) 

Equation (j2.10|) illustrates that for equal mass giant planets, the outer planet's semimajor axis 
must be ~ 30% greater than the inner planet's to ensure that their orbits never cross. The critical 
separation for Jupiter and Saturn gives 02/01 = 1.258, whereas their actual separation 02/0-1 = 1.833 
is much greater and easily satisfies the criterion. 

2.4 The Circular Nonplanar Limit 

In this section we derive a more general form of Eq. ()2.9)1 that incorporates an arbitrarily large 
initial relative inclination between two equal mass planets. The angular momentum components 
(c x , c y , c z ) of a two-body system may be expressed in terms of the angular momentum of the two 
bodies (ci, c 2 ) as: 

c x = c 2x = c 2 sin I sin Q, 

c y = c 2y = c 2 sin/cosfi, (2-H) 
c z = c lz + c 2z = Cx + C 2 cos /, 
where I and Q are the inclination and longitude of ascending node of the outer body. Hence, 

|c| 2 = + /i2# 2 72 + 2/iiAi2^7i72 cos J, (2.12) 
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which differs from the square of the planar expression (J2.7)) only by a factor of cos/. 

The energy of the system is only dependent on the masses and relative semimajor axes of the 
planets. Using the new angular momentum, and setting \i\ = // 2 = //, Eq. (|2.1|) yields a quartic 
equation in 5: 

4 

~ (l + (1 + S 2 + 26 cos/) 2 = 1 + (jj %f - |a* + O (/i 4 / 3 ) . (2.13) 

The right-hand-side is a slowly converging series in mass ratio (the ratio of the first two terms con- 
taining jj, is just [2/i/3] 1//3 ). However, the 1 dominates all terms containing \x for realistic planet-star 
mass ratios, rendering all such terms as minor corrections regardless of their comparable magnitudes. 

Although stability equations such as Eq. (j2.13|) can be solved numerically for 5, an analytical 
solution for arbitrary (not small) initial relative inclination is possible by using computer algebra. 
The exact solution to Eq. (|2.13|) is a series of nested square roots, with radicands of the form 
X + e(I) + /(/, /i), where x is a constant, e is a function of the inclination, and / is a double- valued 
function of the planetary mass and the inclination. En route to deriving Eq. (|2.9)L Gladman (1993) 
encountered radicands of the form x + /(/•*) > an d used binomial expansions based on the assumption 
^ > //. However, making the approximation x ^ H> f° r the radicands with inclination terms is 
insufficient, as e(I) may be comparable to x f° r high inclinations. Further, one cannot make the 
approximation e(I) ^> x; doing so may cause the expansion of the radicand to become indeterminate 
for low inclinations. Instead, we use the following expanded form of the binomial approximation: 

(X + e + /)" = ( X + e) n (l + -J-) n « (X + e) n (l + -^-) . (2.14) 
Repeated application of Eq. (|2.14jl to the exact solution of Eq. ()2.13|1 yields the following result: 



icrit ~ € + 1J\I [ 4+ 1 f e + ^/ i " ) + 



W 3 - 3r/ n\ 



4 + 



cos 2 / 



e + xv^ 3 



+ ■■■ (2.15) 



where 



cos 2 / - cos JV8 + cos 2 /, (2.16) 
cos/ 



V8 + cos 2 / ' 



(2.17) 



X = 3-23-33. (2.18) 

Eqs. ()2.15|) - ()2.18|) provide the critical separation in terms of an initial arbitrary relative incli- 
nation and planetary mass. The equations reduce exactly to Eq. ()2.9)1 in the limit / — > 0°, and 
include a term of the order of unity which is not present in Eq. ()2.9|) . 

Eq. (|2.15J) is in its simplest form because by attempting to simplify it with additional applications 
of the binomial theorem, one would eliminate the /i 1 / 3 term. Further, because the planar limit was 
derived without a small angle approximation, Eq. ()2.15|) allows one to study retrograde planetary 
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motion in addition to prograde planetary motion. If both planets are orbiting in opposite directions, 
their mutual inclination is 180°. Fig. ^illustrates that A crit rises monotonically with inclination. 
Further, the figure illustrates that the Hill stability limit varies little for any mass ratios below 10 -3 . FIG. 1 
In the limit of zero mass and for either perpendicular or retrograde orbits, the Hill stability limit 
reduces to compact rational expressions, as shown in Eq. (|2.19jl . 

lim A crit = 2 + 2V2 lim A crit = 6 + AVs. (2.19) 

i^90° i^l80° 

Although thus far we have considered arbitrary inclinations, hereafter we restrict attention to 
/ < 40°. This bound was chosen based upon the highest inclinations of the known prograde satellites 
in our solar system. Systems exist with prograde satellites that have inclinations relative to the 
ecliptic up to 47°, and are highly unlikely to reside ever in a stable configuration with inclinations 
from 55° < / < 130° (Carruba et al. 2002). 

For most extrasolar planetary systems, the mutual inclinations of multiple planets are unknown. 
However, it is highly unlikely that multiple massive planets could form with very large mutual 
inclinations. Substantial mutual inclination could subsequently develop - for example if one of the 
planets was influenced by a Kozai resonance in a binary (Holman et al. 1997) - but even in such 
cases relative inclinations of tens of degrees appear to be improbable. 

3 Simulation Data 

3.1 Simulation Setup 

3.1.1 The HNbody integrator 

We performed all our numerical simulations with the HNbody integrator (Rauch and Hamilton 
2004 in preparation), which specializes in the integration of systems with a single massive object 
that comprises the vast majority of the system mass. We used the Bulirsch-Stoer integration option 
with the HNbody accuracy parameter set to 10~ 12 (Rauch and Hamilton 2004 in preparation). In 
most cases, energy and angular momentum errors, expressed by (h — ho)/ho and [c z — c Zo )/c Zo , where 
ho is the initial total system energy and c zo is the initial total angular momentum in the direction 
perpendicular to the invariable plane, did not exceed 10~ 7 over the course of the run. In the most 
pathological cases, energy and z-angular momentum errors were conserved to within 10 -4 . Angular 
momenta in the other two directions were typically conserved to two order of magnitudes better 
than the z-angular momentum. 

3.1.2 Initial conditions 

3.1.2.1 Masses, semimajor axes, eccentricities, and inclinations The magnitude of the gravitational 
interaction between two particles results from the masses and relative positions of the particles. To 
facilitate comparison with previous results, we fixed the planetary masses at /ii = /U 2 = 10~ 3 , where 
the subscript "1" refers to the inner planet and "2" the outer planet. Ford et al. (2003) and Ford 
et al. (2001) used these values in their extensive set of simulations, and these values correspond 
well to the Sun- Jupiter mass ratio, a benchmark used in the study of many systems with EGPs. 

Besides masses, the separation is the other factor which determines how planets interact. The 
parameters a, e, and I have the largest effect, while the angular variables have minimal effect. To 
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isolate the effect of inclination, we begin all simulations with circular orbits (ei = e 2 = 0). Further, 
because any two orbital planes in space have just one mutual inclination, we set I\ = 0° without 
any loss of generality. For the equal mass simulations, we fixed initial values of I 2 in multiples 
of 5° ranging from 0° to 35°. We reran the I 2 = 0° simulations with I 2 = 0.001° in order to 
avoid a convergence problem in the code, but will henceforth refer to these simulations as I 2 = 0° 
simulations. Further, because l\ = 0°, the relative initial inclination, denoted by J, will henceforth 
be equal to I 2 . 

The most interesting dynamical behavior occurs for planetary separations that most likely give 
rise to transitory instability followed by quasi-periodic behavior. Based on the definition of Hill 
stability, the outer planet must lie between a 2 = 1 and a 2 = 1 + A crit in order to make orbit 
crossing possible. Therefore, 1 + A crit represented the upper bound for a 2 . The lower bound, the 
boundary at which quasi-periodic behavior largely ceases and at which globally chaotic behavior 
begins, cannot be determined analytically, and no previous empirical estimate has taken into account 
nonzero inclination. For the circular, planar case of two equal masses, Wisdom (1980) found the 
limit (location) of global chaos to be oc /i 2//7 . Not knowing how this estimate holds for nonzero 
inclination, we conducted some preliminary tests, and for ease of interpretation, chose to establish 
an inclination- independent lower bound of a 2 = 1 + A cr i t /2. Having set the range over which a 2 
varies to be 1 + A cr u to 1 + A cr #/2, we chose to decrement a 2 in 300 uniform intervals. 

In Fig. El we verified that for 7 = 0°, all stable systems lie beyond A crit . Therefore, this stability FIG 2 
criterion, although approximate, appears reasonable for the coplanar case. For / > 0°, we found 
that most instances of instability in this separation range correspond to resonances. We explored a 
fixed range of separations (A crit /2 — > A crit ) to facilitate easy comparison between runs. 

3.1.2.2 Angular Variables Consider a planet whose orbit is inclined by / with respect to a reference 
plane. The longitude of ascending nodes establish only the orientation of a system with respect to 
a reference direction, usually the Earth, and do not affect the dynamics of the system. Therefore, 
we set ^1 = ^2 = 0. Having no reason to suspect a priori that particular values of the initial 
mean anomaly would stimulate or inhibit instability, we randomized both M\ and M 2 in the range 
[0, 27r] for each one of the sets of 300 simulations we ran. We repeated each simulation 4 times, each 
corresponding to different (random) mean anomalies. 

3.1.2.3 Time Previous studies (Ford et al. 2001; Adams and Laughlin 2003; Veras and Armitage 
2004) have shown that run times of at least 1 x 10 6 are required in order to allow for dynamical 
settling of the systems. We performed a preliminary set of simulations in order to determine a 
sufficiently long running time for detection of instability, especially in the region of interest, the 
zone outside the global chaos limit. Figure El illustrates the timescales on which systems for four 
different initial inclinations went unstable. The figure illustrates that the number of systems which FIG. 3 
go unstable after 1 x 10 6 trails off appreciably. Specifically, for / = (0°, 5°, 10°, 15°), the percent 

of unstable systems that went unstable within 1 x 10 6 is (91%, 85%, 84%, 74%). This observation, 
combined with CPU time and space constraints, led us to choose 2 x 10 6 as the length of running 
time for all simulations. 
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3.1.3 Classifying instability 

We studied only systems in which both planets remained on bound orbits after 2 x 10 6 ; in 
order to minimize CPU time, we continuously checked for signatures of instability, and when found 
terminated the run. We classified a system as unstable if at any time the semimajor axis or 
eccentricity of either body became negative, or if either body entered the tidal limit of the central 
mass. Weidenenschilling et al. (1984) provide expressions for the tidal limit of two spherical particles 
for a variety of cases. Both tidal shearing and rotation have been neglected in this study. Hence, 
we chose to adopt the largest possible tidal limit, for two spherical bodies in synchronous rotation, 
to our simulations. This condition may be written in our variables as: 

L = 2.29//-5i2 pbnflt , (3.1) 

where L is the tidal limit, and R p i ane t is the radius of the planet, which was taken to be the radius 
of Jupiter. Using Eq. (jH.ljl . we terminated a run when the periapse of either orbit moved within 
the tidal limit. 



3.2 Numerical Results 

We analyze numerical data in three ways: 1) as a whole, 2) in sets of 300 runs, 3) for individual 
runs. We first consider which systems remained stable (defined as both planets remaining on bound 
orbits) over 2 x 10 6 and which did not. Figures 0] and El illustrate the presence or lack of stability 
for all simulations performed. FIG. 4 

Figures E] and El show the stability or instability of each system as a function the initial semimajor ^ 
axis separation and initial relative inclination. As previously mentioned, four trials were conducted 
for each initial set of semimajor axes and relative initial inclination values. Superimposed on the 
plot are the nominal mean motion resonance locations, included to display the correspondence 
between these locations and locations of instability. Note that the horizontal axes are in terms of 
decreasing initial separation, such that the leftmost region of the graph represents the critical Hill 
separation, beyond which the orbits of the planets cannot cross. 

Figures El and El illustrate that systems become unstable either 1) in the global chaos regime or 
2) in the non-stochastic regime at first, second and third-order nominal resonance locations. The 
stable systems of interest are ones where significant dynamical evolution of either or both planets 
has occurred. The greatest measure of this evolution is the maximum semimajor axis deviation 
experienced by a planet. Figures El and [7| isolate systems where either planet's semimajor axis has 
deviated by at least 20% from its initial value, and show that such migrating systems occur either in 
the region of global chaos, or at major mean motion resonances. We chose 20% based on inspection FIG. 6 
of the data; semimajor axis variation for all systems was limited to either a few percent, or easily ^ 
exceeded 20%. We analyze the effective reach of each of these resonances in Section 4.2, but note 
here that each system that undergoes significant radial migration is associated with a resonance. 
Further, in both Figs. El and for the case of I = 20°, unstable systems and migratory stable 
systems associated with the 2:1 resonance tend to lie on the starward side of that resonance. This 
phenomenon could perhaps be explained by the close proximity of the 2:1 resonance to the Hill 
stability limit for / = 20°, as this phenomenon doesn't appear to exist for higher values of I. We 
also point out that Wisdom's (1980) estimation of the global chaos limit, given by the dot-dashed 
vertical lines in Figs. Eland El corresponds well for inclinations up to 15°. 
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We now proceed to analyze individual horizontal rows of data from Figs. 0] and 03 Figure |H1 
illustrates the main features of most of these strips of data by displaying the bounded systems 
(diamonds) which arose from of one set of trials performed with relative initial inclinations of 15° 
(upper panel), 25° (middle panel), and 30° (lower panel). As shown, in most systems the final value FIG. 8 
of the outer planet's semimajor axis differs from its initial value by at most 0.01. In the remaining 
systems, the outer planet may migrate either inward or outward. 

Consider the upper panel of Figure |H1 The system behavior exhibits a clear demarcation at an 
initial outer planet semimajor axis value of about 1.315. This demarcation defines the boundary of 
global chaos. When the outer planet begins life within this boundary, it remains close enough to 
the inner planet that instability occurs and stochastic behavior ensues. When the outer planet lies 
outside this boundary, it usually remains in a quasi-stable configuration, but not always. Consider 
the three conspicuous systems for which a 2 (~ 1.41). Although safely outward of the global chaos 
boundary, these planets show significant radial migration, up to 5. The middle and lower panels 
illustrate how the gross system properties change for greater relative initial inclinations. In the 
middle panel, a few systems migrate at around ~ 1.41, and roughly 10 systems experience migration 
when 1.57 < aiiinitiaV) < 1.59. This same migration regime presents itself in the lower panel, except 
with a large clear gap centered about 1.58. Note in the middle and lower panels, the systems were 
sampled only in the non-stochastic regime. 

The radial locations at which planets migrate are those of mean motion resonances. As the initial 
relative inclination increases, these resonances appear to have greater "reach" , and can more easily 
cause a planet to migrate. This work will focus on identifying the types of resonances influencing 
the systems, but will not delve into the actual interaction at resonance. 

Figure 01 displays the final eccentricity of the sets of systems from the middle panel of Fig. |HJ 
Note that accompanying outward or inward migration is a large and highly variable change in FIG. 9 
eccentricity. Otherwise, the eccentricity remains close to zero, but is nonzero. Such small nonzero 
eccentricities are unreproducible according to Laplace-Lagrange secular theory, as will be shown, 
but may be important in determining the libration widths of eccentric resonances. 

4 Data Analysis 

4.1 Secular Evolution 
4.1.1 Introduction 

Understanding the orbital variations of planets away from resonance enables one to distinguish 
secular dynamical evolution from resonant dynamical evolution. Two planets, initially on circular 
inclined orbits, that settle into quasi-stable orbits around the central mass vary their orbital elements 
according to Lagrange's Planetary Equations, given by Eqs. (14.1)1 for k — 1, 2 (Brouwer and 
Clemence 1961), to a good approximation, 
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dak 2 dTZk 
dt n k a k d\ k ' 
dVt k 1 d7Z k 



dt n k al^l-e 2 k sin 4 dI < 



k 



(4.1a) 
(4.1b) 



dl k -tan|J fe fdK k + dKk\ 1 d7^ ^ ^ 



£ft nfcO^l - e | V 5 A fc n k a\^\ - e\ sin 4 



In Eq. (|4.1|) . 7?.^ represents the disturbing function for either planet, rik the mean motion, X k the 
mean longitude, Wk the longitude of pericenter, and Qk the longitude of ascending node. 

4.1.2 Inclination variation 

4.1.2.1 Laplace-Lagrange theory Two planets on circular inclined orbits that remain undisturbed 
by resonances undergo sinusoidal inclination variations with the same period but different ampli- 
tudes. The goal of this section is to explain this behavior quantitatively, and provide formulas 
for the time variation of the inclination of both bodies in terms of only their masses and initial 
semimajor axes. Laplace-Lagrange secular perturbation theory allows us to derive these analytical 
formulas because of the initial conditions chosen and the approximations used; the results agree 
well with the numerical simulations. 

We apply the secular theory laid out in Murray and Dermott (1999) to our initial conditions 
without resorting to numerical computations. Although this theory relies on several approxima- 
tions, which will be delineated, the results obviate more careful but cumbersome calculations. The 
first approximation sets sin/i « I\ and similarly for the outer planet's inclination. This approxi- 
mation holds for small angles, but even for angles as large as 30°, (I\ — sin/i) / sin/i < 0.05. This 
approximation will be used in the disturbing function and Lagrange's Equations. The second ap- 
proximation is expanding the averaged disturbing function to only second order in the secular terms 
only so that 

TZ secular = - l -abf (If + II) + ab { PhI 2 cos (fix - tt 2 ), (4.2) 

Z 2 2 

Gm 2 . . Gm\ . . , . 

Tvi — \7v secular) , = \7v secular) ■ V±.3J 

a 2 a 2 
Here a = ai/a 2 , the "6" coefficients are Laplace coefficients, and the functional form of the other 
coefficients to second order may be found in Murray and Dermott (1999). The third approximation 
is setting a to its initial value for all time. Inspection of Eq. ()4.1a|) validates this approximation as 
that equation shows that for secular disturbing functions, the semimajor axes remain constant. 
With these approximations, Eqs. (|4.1b|) and ()4.1c|) become, for k = 1, 2, 

dVt k 1 dTZk ,. . N 

—jr = 2T~l^-> 4 - 4a 

dt n k a%l k Oh 

dh 1 m (4.4b) 



dt n k a? k I k dVtj- 
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The inclination vectors pk = Ik sin Qk and = Ik cos allow one to avoid the singularities that 
Eqs. ()4.4j) would cause for low I k and allow one to redefine the averaged disturbing function as, 
using k ^ j, 



(Tl k ) 



n k a 2 k 



-B kk (pi + qt) + B k j (pkPj + qkQj] 



with 



Bkj — Ckj 



— rCtb^P . 

4(m* + m) 2 



(4.5) 



(4.6) 



Here, we have assumed equal mass planets of mass m and a central mass of m*. The values Ckj 
are given in Murray and Dermott (1999). The solutions to this problem, given as inclinations with 
respect to time, are Ik(t) = \fp\ + (jf , where 



p k = ^Tjhj sin (fjt + jj) 

3=1 



(4.7) 
(4.8) 



The quantities 2} are scaling factors, jj are determined by the initial conditions, and fj and Ikj 
represent the eigenvalues and eigenvector components of matrix B, with elements given by Eq. 

4.1.2.2 Application to Inclined Circular Orbits The initial conditions for all our simulations force 
e\ = e% = Q,\ = 0*2 = I\ = 0. In each simulation, we set the relative inclination equal to the 
inclination of the outer planet, which we denote as I . At t — 0, Eq. ()4.7|) then requires 7i = 72 = 0. 
Application of Eq. ()4.8|) . along with our eigenvector choice of In = Ii 2 = 1, gives: 



hi 



J 12 



hll 



22 



I12I2I ' 



T- 2 



Solving for I 2 {f) gives 



h(t) = sjpl + q\ = sjTfll + Till, + 2T{T 2 I 2X h 2 cos (A - f 2 + Tl - 72) t. 
Manipulation of the two eigenvalues fi and f 2 and use of Eq. ()4.9|) shows that 

hit) = J ° - J II + I 2 22 - 2/21/22 cos (B n + B 22 ) t. 
1 + a 2 v 

Insertion of the expressions for the eigenvectors, Ikj-, into Eq. (|4.11|) yields 



h{t) 



1 



1 + y/a 



I Jl + a + 2y/a cos (9t), 



where 



9 = —b 



(i) m -I 



4 2 A/m* + m 



(i/a + a) . 



(4.9) 



(4.10) 



(4.11) 



(4.12) 



(4.13) 
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A similar analysis may be performed to derive I\(t): 



hit) = -^—-I oy /l-cos(et). (4.14) 
1 + v« 

Equations f)4. 12|) - f)4. 14j) exhibit the behavior seen in our secular numerical simulations. The 
negative sign in Eq. ()4.13j) has no apparent physical significance, but does have a mathematical 
significance, and will be needed for later computations. In order to test the robustness of the 
approximation provided by Eqs. (|4.12j) and (|4.14|) . we use a numerical simulation that is most 
likely to give the worst fit to the equations. We obtained such a simulation by performing an 
additional set of 300 runs with a high initial relative inclination (35°) over an initial a 2 range of 
1.54 — 1.63, which is almost centered on the 2:1 (lowest possible order) mean motion resonance at 
1.5874. Therefore, this set of runs samples systems at an initial a 2 increment of 3 x 10~ 4 . The 
lowest <i2 value at which resonantly-induced instability occurs is 1.5544. We use the system with 
ci2 = 1.5541, the closest system to resonance, to test the goodness of our analytical equations to 
secular motions. 

Our numerical results exhibit excellent agreement, as seen in Fig. E3 The inclination range of FIG. 10 
the inner and outer planet, 

< h(t) < -^=, (4.15) 

l + V a 



Io 1 —^ < hit) < /„, (4.16) 
1 + V" 

correspond to those seen in the figure. Note also from Eq. (|4.15jl that because a < 1, the inclina- 
tion of the inner planet, which begins its orbit on the reference plane, exceeds the initial nonzero 
inclination of the outer planet at the maximum value of hit)- Further, Eq. (J4.16)) illustrates that 
the outer planet, which was inclined to the reference plane at t = 0, never exceeds its original 
inclination value. The outer planet oscillates between the reference plane and its initial inclination 
value, but never achieves h = 0°. Conversely, the amplitude of the inner planet is greater, as it 
oscillates between the reference plane and a plane whose inclination value exceeds Iq. 
Eq. ()4.13|) allows one to estimate the period, P, of the oscillations: 



3 



8ny/m + m* . . . 

-af —= ) ±2 m . (4-17) 



mbV ' \V a + a 



All terms in Eq. (J4.17j) except the masses are of order 1, and because m, = 1, the period is 
mostly dependent on the mass of each planet. For Jupiter mass planets, with m = 0.001m*, the 
inclination period of both planets is on the order of 10 3 . 



4.1.3 Node variation 



Having obtained expressions for the inclination of both planets as a function of time, one may 
turn around these expressions and find the longitude of ascending nodes as a function of time. 
Taking f 2 as the nonzero eigenvalue and equating Eq. (|4.7jl with pk = hit) sinQk, then Eqs. 
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Q-fllHJ), and (I ^3D yield: 



sin 9t 



sin 



^2(1 - cos 9t) 



n 2 (t) 



sm 



a sin 9t 



(4.18) 



(4.19) 



/l + a + 2y/a cos 9t 

Comparison with the same numerical simulation that is on the verge of instability, yields Fig. 
im which again shows excellent agreement between the analytical and numerical results. Note FIG. 11 
that the variations in Q are independent of the initial relative inclinations, and that Q\ becomes 
indeterminate during each cycle. However, this indeterminacy is not a result of the approximations 
made in Laplace-Lagrange secular theory. In order to justify this statement, one may use Lagrange's 
equations in their full generality to obtain an alternate expression that also exhibits this degeneracy. 

The gravitational interaction between two planets, when not in a strong resonance, can be 
considered by disturbing functions that incorporate only secular terms. Hence, with the inclusion 
of secular terms up to second order in inclination, we use: 



{Tlx) = [n x ) secular = C + d (si + s 2 2 ) + C 2Sl s 2 cos (fi 2 - fix), (4.20) 

where s\ = sin (ii/2), s 2 = sin (I 2 /2), and Co, C±, and C 2 are constants. Application of Eq. ()4.2()|) 
to Eq. ()4.1b|) (setting e\ — without loss of generality) yields, with £ = Q 2 — fi 1? 

dQi C\ C 2 s 2 cos£ \Ci\ 
dt 2n\a\ 4niaf S\ 2ri\a\ 

The last equality follows from the relation between C\ and C 2 (Murray and Dermott 1999). The 
expression for ^2 is equivalent, except for an interchange of subscripts. Eq. (14.21)1 illustrates that 
as Ii approaches zero, Cli becomes indeterminate. Further, Eq. (j4.21j) illustrates that the condition 
for the longitude of nodes to increase may be expressed as: 

Cl% > when 
Q 2 > when 
Further, differentiating Eq. (J4.18)) gives 

(hit) = --, (4.23) 

a constant. Eqs. ()4.22)1 and ()4.23jl provide a mathematical basis for the time rate of change of VL\ 
seen in Fig. ^2 Differentiation of Eq. (|4.18|) yields a somewhat, less cleaner result (than Eq. 14.23)1 
that contains a to take into account the curvature seen in Fig. EI Further, integrating Eq. ()4.23)) 
directly illustrates that Q± is a linear function of time, but leaves an undetermined constant that is 
absent from the exact Eq. (|4.18)1 . 



s 2 cos £ 
si 



(4.21) 



■S2COS4 

^ t, 

si 
Si cos£ 



> 1. 



S-2 



(4.22) 
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4.1.4 Eccentricity variation 

Unlike for a planar eccentric system, in which 1 = for all time, for a nonplanar circular system, 
e 7^ for all time. Hence, both planets develop a nonzero eccentricity, which can play a crucial role 
in determining the type of resonance acting in a system. However, second-order Laplace-Lagrange 
secular theory cannot reproduce these nonzero eccentricities, instead predicting that e = for all 
time. This result, which may be presupposed based on the decoupling of eccentricity and inclination 
in second-order Laplace-Lagrange theory, is often made heuristically using degeneracy considera- 
tions. In this section, we confirm this result rigorously, linking the eccentricity and inclination 
through the Lagrange equation for dfl/dt. Having developed the formalism to derive analytical 
formulas for the circular inclined case, we then briefly show how to obtain similar equations for the 
coplanar eccentric case, and the complications that result from doing so. 

4.1.4.1 Forever circular orbits The Lagrange Equation given by Eq. (l4~Tb|) is an exact formula, 
derived without approximation. We will use this formula to solve for the eccentricity as a function 
of time, as we now know all the other quantities in the formula from secular theory. These "other 
quantities" are in good agreement with the numerical results. The secular part of the disturbing 
function expanded to second order in both eccentricities and inclinations may be written as (Murray 
and Dermott 1999): 

^secular = ^ / 2 + £ T\ _ * & (1) / fi 2 + g\ + ^ ^ _ ^ + ^(D^ cog _ ^) ( 4 . 24 ) 

Z 2 2 

where K\ and K 2 are constants that will be irrelevant to the subsequent analysis. Equation (|4.24|) 
may be used in the expressions for both disturbing functions of the problem, such that (Murray 
and Dermott 1999): 

K x = — a7^ cular K 2 = ^i a ^ cular . (4.25) 
Taking the partial derivative appropriate to Eq. (j4.1b|) gives: 

dlZl Gm 2 a%% /l V 

cos -ii |s 2 cos4 — si\ . (4.26) 



dh 2a x \2 

Inserting Eq. ()4.26|) into Eq. (j4.1b|) and solving for e\ gives 



e x (t) 



Gm 2 a 2 b^ 2 [s 2 (t) cos^(t) — s 1 (t)\ 
\ \ 4afVm* + mti 1 (t)s 1 (t) 



(4.27) 



Similarly, 



e 2 (t) 



(4.28) 



Gmiab^ 2 [s\(t) cos^(t) — s 2 (t)] 

\ \ 4a| A/m* + mtl 2 (t)s 2 (t) 

Equations ()4.27|) and ()4.28|) illustrate how both planets' eccentricities vary through time due to 
secular effects alone. Knowing an analytic expression for the eccentricities at a particular time, one 
may now proceed with a secular analysis similar to the one performed for the planets' inclination. 
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For a disturbing function that includes eccentricity terms Ajk, and for scaled eccentricities h k = 
efcSinti7 fc and K k = e k cos uJk, one may write down a set of equations similar to Eqs. (|4.6|) - ([4.8|) 
(Murray and Dermott 1999), 




(4.29) 




(4.30) 



hk = ^ Sjtkj sin (gjt + fy) 

3=1 



(4.31) 



2 



Kk = 



Sje kj cos (gjt + Pj). 



(4.32) 



Note that because of the lack of "mutual eccentricity" (like mutual inclination), the coefficients 
Ajk lack the symmetry of the Bj k terms, which all had the same Laplace "6" constants. This pre- 
vents significant simplification of the roots in the eigenvalue {gu) and eigenvector (e^) expressions. 
Further, the phase angles j3 now are nonzero because of the temporal conditions used. 

Earlier we could not apply Laplace-Lagrange secular theory to the eccentricities, because initially 
ei = e 2 — 0, and hence the secular equations were degenerate at the only known temporal condition. 
Now, however, with Eqs. ()4.27|) and ([4.28)1 . we may obtain another time condition (not necessarily 
an initial condition). Denote by to the time at which the inclinations of both planets are equal, a 
condition which occurs twice in every inclination cycle. Equations (|4.12j) and ()4.14|) then show that 
to satisfies 



Evaluating ei(t ) using Eqs. (j4~T2l - (l4~T4l . lOSt - ljCTty . (jPSJ) . (jQ7jl . and (g~35j) yields, after 
much algebra, 



One may also show that e 2 (to) = 0- Hence, because the eccentricities are zero at both t = 
and t = t , one may use Eqs. (|4.31|) and ()4.32|) to illustrate that /3i = /3 2 , an d hence g\ = g 2 . 
However, as shown by Eqs. ()4.29|) - ([4.30p . the two eigenvalues of matrix A do not satisfy g\ — g 2 - 
This contradiction forces both eccentricities to be zero for all time, an unphysical result. Hence, 
secular theory cannot predict the eccentricity evolution of two planets on initially circular orbits, 
even if those planets do not ever drift into a strong mean motion resonance. 

4.1.4.2 Planar eccentric case Although not directly relevant to the results of this work, an ex- 
pansion of the previous results to the planar eccentric case could perhaps prove useful for future 
studies. Just as previously we set the eccentricities and nodes to zero, here we set Ji = 1% — vd\ — 
xz>2 = 0. However, unlike inclination, eccentricity is not defined to be a relative quantity. If we set 
ei(t = 0) = and e%{t = 0) = e2o, then 




(4.33) 




(4.34) 



S 



i — 



ene 22 - e 12 e 2 i 



e 2 oei2 



S 2 




(4.35) 
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Otherwise, if we set e 2 (t = 0) = and ei(t = 0) = e w , then 

eioe 2 2 



Si 



ene 22 - e 12 e 2 i 



S 2 = -S^. 

e 22 



(4.36) 



Solving for e k {t) gives 



e k (t) = + n 2 k = ^Sfe 2 kl + Sle\ 2 + 2S 1 S 2 e kl e k2 cos (gi -g 2 + Pi~ &) t. (4.37) 

Regardless of the choice of a nonzero e w or e 20 , Pi = P 2 from Eq. (|4.31|) . Hence, one may 
obtain expressions with the same functional form as Eqs. (j4.12|) and (j4.14|) for eccentricity. Direct 
computation of the eigenvectors gi and g 2 (Murray and Dermott 1999, p. 318) shows that 



9i ~ #2 



2 



\ 



( m2 V _ 2— 

\rrii J nil 



b) 



(2) 



3/2 

w 

J 3/2, 



a? + a 



(4.38) 



The asymmetry produced by computing eccentricity instead of inclination is contained entirely 
within the ratio of Laplace coefficients in Eq. (|4.38|) : when this term becomes unity, and in the 
limit of equal mass planets, Eq. (j4.38J) reduces to Eq. ()4.13|) exactly. Thus, the deviation in the 
period of oscillations in the inclined circular case to that from the planar eccentric case goes as, 
when a<l: 



1 - 



l(2) 
2/2 

h/2 



1 — -a 
4 



(4.39) 



where polynomial expansions for the Laplace coefficients in a can be found in Murray and Dermott 
(1999). 



4.2 Libration Widths 

By inspection of Figs. |U|7l one might conclude that all unstable or significantly evolving systems 
that initially lay outside of the global region of chaos are associated with a resonance. This section 
will attempt to support this hypothesis by classifying the likely resonances that are interacting, and 
quantifying the effective "reach" of such resonances. 

Two bodies are said to be "captured" into a resonance when one or both of the bodies drifts 
closer to the other until a commensurability is attained. If the planets are moving away from each 
other, resonant capture cannot occur (Sinclair 1972; Yu and Tremaine 2001). However, resonant 
capture is not assured if two planets approach each other. The perturbations on two planets in 
quasi-stable orbits alone does produce instances where the planets are momentarily drifting toward 
one another, but for our simulations, this drift was not sufficient to cause any resonant capture. 

Just because a planet's initial semimajor axis lies close to a strong resonance does not mean the 
planet will achieve a state of near-exact resonance with that particular commensurability. Further, a 
planet may undergo a series of resonant encounters, each on a different timescale, so as to obliterate 
any record of its previous locations. These considerations should be remembered when comparing 
the initial and final semimajor axes of the planets studied here. Further, a planet has zero probability 
of achieving a particular exact resonance, but a finite probability of drifting inside a relevant interval 
centered on an exact resonance. This interval is most typically classified as the "libration width", 
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and is relevant in the sense that a planet drifting inside may suffer a significant radial excursion if 
not captured by the resonance. 

The gaps in the strings of black diamonds in Fig. |S] appear to increase with the initial relative 
inclination. These gaps visually may be thought of as libration widths. In order to affirm and 
quantify this trend of increasing libration gap width with inclination, we now proceed to derive an 
analytical estimate of the gap width. Murray and Dermott (1999) provide such an estimate for the 
eccentric, planar case with one negligible mass; we consider the nonplanar, non-circular case with 
no assumptions about the masses. 

Suppose one wishes to find the libration width of a resonance whose resonant argument is <fi = 
ji Ai + j 2 \ 2 + J3^i +34&2 + J6^2- One of the disturbing functions may then be expressed as: 

= 9H^ [{ni)secular + C^^s^s? cos0]. (4.40) 
a 2 

A complete description of the nontrivial computation of C\ may be found in Murray and Dermott 
(1999). Application of Eq. (ETTal) yields 

^ = — = - 2 ^^C 1 et l ef^s^ sin0. (4.41) 

at n\d\ 0X1 n\a\a 2 

which may be rewritten as 

h x = 3 ^ 2il gi4 i3 'e2 4 '4 i5l 4 i61 sin0. (4.42) 
afa 2 

Similarly, the other disturbing function for the problem can be used to write: 

n 2 = ^!^ C 2 e' J3| e' 34l 4 i5l 4 J61 sm0. (4.43) 
a 6 2 

Ignoring the second derivatives of angular variables, one obtains: 

" hn 1+ hh 2 = ef^sf^ (^P^C, + 3 -^C 2 ) sin0. (4.44) 



a 2 a 



2 



Equation ()4.44j) is the equation of a pendulum, where the angular frequency may be repre- 
sented as the coefficient of sin </>. Comparing the maximum to minimum energy of this pendulum, 
and relating the result to Eq. ()4.43|) yields: 

dn 2 = ±— sin ( -(f)) dtp. (4.45) 



32 V2 

Integrating Eq. ()4.45|) and converting to a 2 gives the libration width, 5a 2y iib. 



, n ,4o2u;o 4oa / b3 | j J4 | b5 | j J6 | [ SGm.jf 3Gm x3 \ \ \M M M \M 

da 2 ub — ±- — — — ±7- — -4/ei e 2 s l s 2 — 5 H o — L 2 ) on e 1 e 2 s x s 2 ■ 

3n 2 j 2 3n 2 j 2 ]j \ afa 2 J 

(4.46) 

Equation ()4.46|) compares favorably with a similar expression derived by Champenois and Vienne 
(1999b) to study the Mimas- Tethys 4:2 resonance. One observes that 5a 2 jib monotonically increases 
with an increase in the eccentricity or inclination of either body, and for two equally massive 
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planets is ~ \[2 greater than that of the restricted three-body case. Further, the d'Alembert 
relations require that the order of a resonance equals |j 3 + j'4 + j'5 + and that the quantity 
|j5+j6| be even. Hence, any first and third-order resonances must include a contribution from 
one of the planetary eccentricities. Therefore, the inability of secular theory to reproduce the 
observed eccentricity evolution unfortunately shrouds one's ability to determine which resonances 
are acting. Further, because the eccentricities can reach in excess of 0.01, their contribution to 
the libration width often rivals the contribution of the inclinations, preventing one from neglecting 
the eccentricities in resonance calculations. For second-order and higher-order resonances, with 
planetary eccentricities of 0.03, an initial relative inclination of 10° produces libration widths that 
are comparable to those produced by eccentricities. Any higher initial inclination will cause vertical 
resonances to dominate. 

Being functions of eccentricity and inclination, libration widths should also be a function of time. 
In Fig. ^1 we provide a sample of the magnitudes of the libration widths for two independent-of- 
eccentricity 5:3 vertical resonances (the "mixed" resonance with j'5 = j§ = —1, and the "unmixed" 
resonance with either j 5 = —2 or j 6 = —2), using Eqs. ()4.12|) and (j4.14|) to model the time evolution. 
As shown in Fig. El an unmixed vertical resonance is most likely to take effect in the middle of FIG. 12 
a dynamical inclination period, and a mixed vertical resonance is most likely to take effect at the 
beginning (or end) of a dynamical inclination period. 

The large phase space covered by our simulations necessitated an undersampling of the data, 
preventing us from exploring the detailed interactions at resonance. Figures EE illustrate that the 
2:1 resonance location affects a relatively large number of systems. The 2:1 resonance is a first- 
order resonance, implying that planetary eccentricity, only, determines the dynamical evolution of 
systems residing in that commensurability. However, the 2nd-order 4:2 resonance, 3rd-order 6:3 
resonance, or such higher order resonances may be acting instead, and thus inclination may play 
a role as well. Only detailed, small-timestep calculations of custom-models of such resonances can 
further pinpoint the dynamics induced there. 

Table 1 documents the effect that every first through third-order resonance which overlapped with 
the phase space studied has on all systems that exhibited significant dynamical evolution (where 
either semimajor axis varied by over 20%). Of the 9600 systems modeled, 87 stand out because TABLE 1 
they: 1) remained stable over 2 Myr, 2) began outside the region of global chaos, 3) and have shown 
significant radial migration of at least one planet. Table 1 illustrates that 85 of these 87 systems 
(98%) had initial configurations where both planets lay within two libration widths of a 1st, 2nd, 
or 3rd-order resonance. In the table, columns labeled ±lu> and ±2w list the number of systems 
that initially lay within one and two libration widths respectively of the resonance location in 
question. The libration widths used were the maximum widths (computed from Eq. I4.46j) attained 
over the secular evolution of the systems, which were modeled using Eqs. (14.12}) and ()4.14|) for the 
inclinations and with fixed orbital eccentricities of 0.03 (an upper bound on the eccentricities seen in 
the simulations during secular evolution) for each planet. The rightmost column displays the percent 
of systems exhibiting significant radial migration that initially lay within two libration widths of 
any lst-order thru 3rd-order resonance for given relative initial inclinations, /. For I = 20° and 
/ = 25°, some systems resided within twice the libration widths of both the 8:5 and 5:3 resonances, 
or both the 5:3 and 7:4 resonances. The bottom row shows that 84% of systems that fulfilled the 
three criteria above were affected by the 2:1 or 5:3 resonances. The two systems not "picked out" 
by the conditions stated above initially lay within 0.01 of the fourth-order 7:3 resonance, located 
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at 1.759. This location represents the largest outer planet semimajor axis for which two planets on 
inclined, circular orbits exhibit radial migration, as Table 1 illustrates that despite residing within 
the region that allows crossing orbits, the stronger but further out 5:2 and 3:1 resonances do not 
induce such behavior. 

5 Mapping Results 

5.1 Background 

We now investigate the extent to which encounter maps can be manipulated to reproduce the 
behavior seen in the integrations. Encounter maps represent useful tools for describing the geometry 
of successive conjunctions between a test particle and a massive perturber, one which is orbiting a 
more massive central body. These maps also have been applied to the orbits of planetesimal swarms 
(Namouni et al. 1996, henceforth referred to as NLTP) and outer solar system objects (Duncan 
et al. 1989). The primary advantage of maps over numerical integrations is the 10 to 100-fold 
decrease in computation time. Although these maps traditionally model systems with a massless 
test particle, we show that they may be used to model systems with two equal mass planets, and 
compare the results of these maps with our numerical integration results. 

Duncan et. al (1989) derived an encounter map for the case where the inclinations of the test 
particle orbit and perturber orbit are zero, the initial eccentricity of the test particle is small, and the 
initial eccentricity of the perturber is carried out to first order. Their first-order map used the small 
eccentricity approximation given by Julian and Toomre (1966). NLTP derives a more general map 
that allows for initial eccentricities and inclinations (which will be henceforth denoted by lowercase 
"i", in order to maintain consistency with NLTP's variables) of the perturber to be expanded to 
arbitrary order about e = and i = 0. Despite this expansion about zero inclination, we show 
that NLTP's map, to second-order only, can faithfully reproduce the qualitative gross features of 
the relatively high-inclination systems we have studied. With this map, we then extrapolate to 
semimajor axis regions of interest that we did not explore with the HNbody code due to CPU 
constraints. 

5.2 NLTP's Map Generating Algorithm 

NLTP derives a planar, second-order encounter map, by using the variables and equations we 
summarize in Table El In order to maintain consistency with the rest of this study, we define TABLE 2 
the quantities they use as the mass m, semimajor axis a, eccentricity e, inclination i, longitude of 
perihelion u, and longitude of ascending node Q. They further define the mean longitude A, true 
longitude 8, Hill parameter e, and action angle variables J, J, and K. The Hill parameter scales a, 
e, and i into Hill units, which enable NLTP to derive the desired mapping algorithm. 

The quantity NLTP use to perform expansions in eccentricity and inclination out to arbitrary 
order is an effective potential, defined by W . This potential is derived from Hill's equations, where 
Hill's equations are obtained from the following approximations: 

TO2,m 3 <TOi, e 2 ,e 3 <l, 
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W takes the form of a Fourier series, where, with j = y — 1, n G Z, 

oo 

W = Wne^-^, (5.1) 



O /-7T 

_ / g-in(|e* sin(t-£j)+t-£j) 



<27T 



(2n i\ 
y ((e* cos(t - u) + l) 2 + z 2 cos 2 (t - n)) 5 J d(t 



(5.2) 



Here i^o is the zeroth-order modified Bessel function of the second kind. Equation (|5.2j) gives 
the exact form of W n . For details involved in the derivation of this "interaction potential", we refer 
the interested reader to NLTP. 2 In order to use different orders of approximation, NLTP expand 
this equation in the following form: 

W n = \ J2 W™<*i? + ... (5.3) 

\n\<p+2q<M 
parity (n) =parity(p) 

The quantity W%' q is expressed in terms of Bessel functions and the difference, (fl — u). The 
terms W% ,q are generally a function of (fl — uj) only, but are constants when q = 0. For < q < 4, 
WP' q is a function of Q only. Values for many of these terms are given in Appendix B of NLTP, 
and are not be reproduced here. M represents the order of the approximation of the desired map, 
while p and q are integers that do not exceed M. Hence, the second-order planar map that NLTP 
derive corresponds to M = 2, q = 0. 

The effective potential W is used to derive the encounter map from the following equations: 

d 

In+i — In + -7rrW(uj n , fl n , X n , I n +1, Jn+1, Kn+l), (5-4) 
0(p 

d 

Jn+1 = Jn + -7r-W(uJ n , Q n , A n , I n +1, Jn+1, Kn+l), (5-5) 

dip 

Kn+l — In+1 + Jn+1 ~ In ~ Jn + K n , (5-6) 
8 

^n+1 = ^n — 777^(^1 ^n, ^n, In+1, Jn+1, Kn+l), (5-7) 

ol 

d 

fln+l — ^n — 777^(^71; ^rn A n , / n +l) Jn+1> ^n+l)j (5-8) 



A n +1 — — ^T~7W(LU n , f2 n , A n , J n +l, J n +1) Kn+l) + TT \ cTT? ' (5-9) 

<9il 3e y 8il n+ i 

Note that I n +i, J n +i an d K n+ \ appear implicitly in the above equations. As a result, solving 
for the complete map analytically in closed form is impossible for all but the lowest order approx- 
imations. Therefore, NLTP "keep only the free motion contribution to the longitude" by setting 

2 One should be aware of two typos in NLTP: the cosine and sine in their Eqs. (37) and (38) should be interchanged. 
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a r = a>2 in Eqs. (|5.4|) - (j5.9j) . Additionally, they dispose of the differential term in Eq. (|5.9|) be- 
cause they claim that this approximation holds when the bodies, "do not suffer large excursions in 
space." Because we are using NLTP's algorithm only to aid in our effort of quantifying what initial 
conditions produce instabilities, and not to explore the evolution of any one particular system, we 
zero out the differential term in Eq. (j5.9|) as well. 

One can see that from Eq. (|5.H|) that when W n is expanded, many terms with a linear dependence 
on eccentricity will result. These terms, from the definition of the action-angle variable J, are 
proportional to y/1. Then from Eq. (|5.7j) of the map, one sees that the derivative of W is taken 
with respect to /, resulting in singularities for the aforementioned terms. In order to circumvent 
these singular terms, NLTP first splits the potential, W, identifying the portion that contains the 
singular terms as W S i ng . NLTP then iterates W S i ng with the symmetry-preserving variables, h and 
k, to create an intermediate step between n and n + 1 that removes any degeneracy. With primed 
variables representing the intermediate step, the resulting equations read: 



h' = h n + —W smg (k n , n n , A n , ti, J\ K'), (5.10) 
d 

k — k n -Qj^W S i n g(k n , Q n , X n , h , J , K ), (5.11) 

, h' + k' , 

r = — ^— , (5.i2) 

9 

J = Jn + -q^W s i n g{k n ^ Q n , A n , h , J , K ), (5.13) 

K> = +k 2 K K + J' - J n + K n , (5.14) 

to = arctan — , (5.15) 

n' = n n - ^w smg (k n , n n , A n , ti, f, k'), (5.16) 
d 

A' = A n — ~^zpWai n g(k n , Q, n , A n , h , J , K ). (5-17) 



Similarly to Eqs. (j5.4j) - (j5.6|) . Eqs. (|5.10j) - (|5.17jl are implicit functions, and thus inhibit one's 
ability to write down explicit analytic expressions for higher order maps. Equations ()5.4|) - (J5.17)) 
represent the algorithm one may use to derive NLTP's encounter map. The number of terms kept 
in the Taylor expansion of W (Eq. 15. Hj) determines the accuracy of the map. We now use this 
algorithm to derive a second-order spatial map that is applicable to the systems integrated in this 
work. 

5.3 The Nonplanar Second-Order Map 

NLTP use planar maps of up to tenth order, and provide the full expression for the second-order 
eccentric planar map, but do not study nonplanar maps. In this section, we use NLTP's general 
mapping algorithm to create a second-order eccentric nonplanar map, which correctly reduces to 
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NLTP's second-order eccentric planar map in the limiting case of zero inclination. Then we show 
that with a modification, the second-order eccentric nonplanar encounter map roughly reproduces 
the same areas of instability seen in the numerical integrations. We finally apply this map to areas 
of semimajor axis phase space that remained unexplored by the numerical integrations. 

The six degrees of freedom used to initially define planets for these maps is the set {a, e, i, A, u, fl}. 
This is the same set of variables we used for numerical integrations, except here we use the mean 
longitude (A) instead of the mean anomaly (M). Given initial values of these six variables, we can 
define the moments Iq = eg/2, J = 2q/2, Ko = 3<2q/8, h = y/2I cosc<j , and k Q = y/2I sin u> , then 
iterate with our second-order nonplanar eccentric map, which we obtain after extensive manipulation 
of Eqs. (|5.4|) - (j5.17|) . The final map reads: 



ti = h n + 



2W£'° sin A n 
ai 



<3 



2W 2 1,0 cos A n 



k — k n 



n o 

a 3 



_ ti 2 + k' 2 

2 ' 
k' 

uj = arctan — . 

h' 



j Jn 



5.18) 
5.19) 
5.20) 
5.21) 

5.22) 

5.23) 
5.24) 

5.25) 
5.26) 

2^ ,j. - ,, - \5.27) 

Some features of Eqs. ()5.18|) - (j5.27|) deserve discussion. Here, 1/2 represents i^i(4/3), where Ki is 
the modified first order Bessel function of the second kind. The sine term in Eq. (15.23)) is a compact 
way of expressing the four double angle terms that result from the application of Eq. ()5.4|) . Recall 
that the W% ,q in the mapping equations are constants resulting from the Taylor expansion of W 
about e* = and = 0. Note that Eqs. ()5.18j) - (j5.27j) contain only the coefficients W 2 '° and W^ 2,0 , 
which represent coefficients of zero inclination terms (see Eq. 15. 3|) . The coefficients of the inclination 
terms in the expansion do not vanish, but rather were combined into trigonometric arguments of 
angle differences, which are manifested in Eqs. ()5.22|) . ()5.23|) . and ()5.27|) . 



1 + m a ~ *y 2 cos[2(A n - u')] sin[2(fi n - lo>)\ 

_ V_ - fJ n+1 y 2 af sin (2^ n - A<j>' + 2A W ) 
n+1 ~ 1 + 8 sin [2(A n - w>)\ afW*' 



... /sA".. , , \ ~ 
A„ +1 = A n + -^ 3 J , 

Un+1 =J- [4W 2 2 'V^ 3 cos [2(A n - uo')} - 2a^ 3 ] 



-y 2 cos[2(fi n - (v'))a 3 3 cos [2(A n - uj')} + 2a 3 3 
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Applying the map given by Eqs. ()5. 18)) - ([5.27)1 to two equal mass planets provides a null result, 
unless the Hill parameter is modified suitably. Because the equal mass problem deals with two 
bodies whose orbital parameters change at comparable and non-negligible rates, the parameters 
in the map must be scaled accordingly. Doubling the Hill parameter satisfies this scaling, so that 
e = 2 [(mi + m 2 )/m ] 1 ^ 3 . Applying this new Hill parameter to the second-order nonplanar map 
yields Fig. El which can directly be compared with the instability seen in the numerical integrations 
for all runs with I = 0° — 15° with the mapping results. The map mirrors the numerical integrations FIG. 13 
to a rough extent, providing an estimate for the limit of global chaos and for instability due to 
resonance. The map was run with initial, random mean longitudes, and was run for just « 1000 
conjunctions. Also, five times as many initial semimajor axis values were sampled, but otherwise 
all initial orbital parameters were the same as the numerical integrations. Unstable systems were 
flagged as those systems whose relative eccentricity exceeded 1. Such systems presented themselves 
in much less time (both real time and CPU time) in NLTP's encounter map than in the actual 
numerical integrations. Although NLTP's map provides a fast, convenient way to discover instability 
in a system, we find that the results are not accurate enough with respect to our integrations to 
warrant an analysis of individual runs. Hence, we use the map only to explore the gross properties 
of systems. 

In Fig. El the map output reproduces the instability arising from the 5:3 resonance. However, 
the instability is offset from the actual location of the resonance. We assume this offset results from 
the Taylor expansion of the map around i = 0; as the planet's initial relative inclination increases, 
the map is less able to reproduce the actual behavior of the systems. Applying the map to relative 
initial inclinations of 25° and higher results illustrates that the map fails to reproduce the instability 
at the resonances seen in Fig. However, for inclinations < 15°, the map can provide a qualitative 
estimate of the incidence of stability. Given this conclusion, we now proceed to consider regions of 
phase space unexplored by the numerical simulations. 

Having established that NLTP's map can produce qualitative features of the instability seen in the 
integrations, we ran the map for a wider range of phase space than that covered by the integrations. 
Figure ITU shows the results for a sampling of 6, 000 equally spaced outer planet locations over the 
interval [1, 2] for 6 different inclinations. One observes that most systems outward of the limit FIG. 14 
of global chaos, and the vast majority of systems outward of the Hill stability limit, display no 
significant dynamical evolution, with the notable exceptions of the 5:3 resonance systems and a 
hint of instability at the 2:1 resonance location. An equivalent exploration of phase space for the 
interval [2, 3] was performed. Not a single one of those 36, 000 runs became unstable. A final 
exploration of [1, 2] phase space was performed for inclinations higher than 15°, in order to test the 
robustness of the Taylor expansion about % = 0. The results are presented in Fig. El Inclinations FIG. 15 
greater than 25° illustrate nothing discernible, as they fail to reproduce the stable gap located at 
1.35 produced by lower inclination trials. 

Therefore, encounter maps 1) successfully reproduce the gross properties of the instabilities seen 
in the numerical integrations, suggesting that such maps can be used to model two massive planets, 
and 2) illustrate that no instability was seen when a 2 > 2, lending support to the concept of Hill 
stability. 
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6 Discussion 



For application to known extrasolar planetary systems, the first question is whether multiple 
giant planets are typically formed in essentially coplanar configurations, or instead have significant 
mutual inclinations. Observationally, this is not known, although in some known multiple systems 
stability considerations constrain the maximum mutual inclinations of the planets. For example, the 
two planets in 47 Ursae Majoris cannot lie on orbits inclined by more than about 40° to each other 
(Laughlin et al. 2002), while in the more complex three-planet system around upsilon Andromeda 
high inclinations of the outer planets appear to be likewise disfavored (Lissauer et al. 2001; Rivera 
2001). These limits, however, while ruling out extreme scenarios such as counter rotating orbits, do 
not preclude mutual inclinations large enough to significantly alter the subsequent orbital dynamics. 
Moreover, direct observations of dust debris disks, especially that around f3 Pictoris, reveal warps 
of the order of 10° which extend inward to radii where giant planets orbit in the Solar System 
(Wahhaj et al. 2003). Although in these cases the observed warps may be the consequence of 
planet formation rather than a relic of the initial conditions, there are also indications that warps 
might be present — at least at large disk radii — at earlier epochs when the protoplanetary disk was 
still gas-rich (Terquem and Bertout 1993; Launhardt and Sargent 2001) 

If gravitational scattering is important in the early evolution of extrasolar planetary systems 
(either in initially coplanar configurations, or in the inclined systems considered here), then it will 
influence both the orbital radii and eccentricity of the observed post-scattering planets. Inward 
migration as a consequence of two planet scattering is unlikely to explain the existence of hot 
Jupiters at orbital radii a < 0.1 AU, since the progenitor multiple systems would themselves need 
to have formed at small radii where massive planet formation is probably inefficient (Bodenheimer 
et al. 2000). A more promising application is outward migration from radii comparable to that of 
Jupiter to several tens of AU, where standard core accretion models (Pollack et al. 1996; though 
see also Goldreich et al. 2004) predict formation timescales that are comparable to or longer than 
typical protoplanetary disk lifetimes. There is plausible circumstantial evidence that a population 
of such planets exist, from observations of asymmetries in several debris disks (Wilner et al. 2002; 
Quillen and Thorndike 2002; Holland et al. 2003). If planet-planet scattering is common in the 
early evolution of planetary systems — as would be required if this is the origin of the typically 
substantial eccentricities of extrasolar planets — then the results presented here suggest that the 
large radii planets could be a byproduct of the same process. Large eccentricities of the outward 
scattered planets would then be common, a prediction which appears to be consistent with the 
limited inferences possible from existing data. 

In addition to the (small) fraction of known extrasolar planets that have been found at sur- 
prisingly small orbital radii, a large fraction of the planets further out have significantly eccentric 
orbits. For massive planets outside the hot Jupiter zone, eccentricity is distributed approximately 
uniformly in the range < e < 0.7 (Marcy et al. 2003), with a smaller number of outliers at even 
higher eccentricity. This is certainly qualitatively consistent with a high frequency of planet-planet 
scattering, though whether multiple planet formation is required to yield significantly non-circular 
orbits depends on the rather uncertain issue of whether a single planet's eccentricity can be excited 
as a consequence of interaction with the protoplanetary disk (Artymowicz 1992; Papaloizou et al. 
2001; Goldreich and Sari 2003; Ogilvie and Lubow 2003). Encouragingly for the multiple planet 
scenarios (Rasio and Ford 1996; Ford et al. 2001; Marzari and Weidenschilling 2002), two-planet 
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integrations that include a realistic range of planet masses (Ford et al. 2003) provide a better 
quantitative match to the observed eccentricity distribution than did prior simulations assuming 
equal mass planets. 

Overall then, the available evidence is consistent with a high frequency of gravitational scattering 
in initially multiple massive planet systems. The calculations presented here, however, confirm that 
gravity is unlikely to be the only significant agent at work. If planets are formed with extremely small 
orbital separations, then interesting dynamical behavior (excitation of eccentricity and/or ejection 
of one of the planets) is inevitable, but the timescale is typically short. A significant fraction of 
these eventually unstable systems develop instability in 10 2 — 10 5 years, a shorter timescale than 
the typical disk dissipation timescale (Hartigan et al. 1990; Wolk and Walter 1996). This implies 
that purely gravitational scattering experiments are physically inconsistent — in reality either such 
unstable initial conditions could not arise, or additional forces such as torques from the remnant 
protoplanetary disk would have to play a role. If, alternatively, multiple planets are formed with 
larger separations that exceed the global chaos limit, then as we have seen the most interesting 
dynamical behavior is restricted to the vicinity of low-order mean motion resonances. If planets 
form at random locations with the protoplanetary disk, then only a few percent of all systems 
would satisfy this condition. In fact, approximately ten percent of known extrasolar planets are 
found within multiple systems, and these planets have a similar distribution of eccentricities to 
those planets with no (known) companions (Marcy et al. 2003). Several of the multiple systems are 
either definitely or probably in resonance. Taken as a whole, we believe that these properties are 
consistent with a scenario in which both disk-driven migration and planet-planet interactions play 
critical roles in explaining the origin of extrasolar planet eccentricities (e.g. Snellgrove et al. 2001; 
Lee and Peale 2003; Chiang 2003), although quantifying the statistical properties of the resulting 
planetary systems within such models remains difficult. 

7 Conclusions 

We have explored the dynamics of two close massive planets on inclined, nearly circular orbits 
using both analytical and numerical techniques. Our findings are summarized as follows: First, our 
main result is that significant outward radial migration of EGPs through gravitational scattering 
alone is possible and occurs either inside the region of global chaos or at the first, second, or third 
mean motion resonances (especially the 2:1 and 5:3 resonances). The few cases of migration that 
do not occur at such resonances migrate at the next lowest-order resonance locations. Scattering 
can also lead to some inward migration, but for two equal masses, the degree of inward migration 
is constrained by energetic arguments to be modest (Lin and Ida 1999; Papaloizou and Terquem 
2001; Adams and Laughlin 2003). 

Our ancillary results are 1) the Hill stability limit for arbitrarily inclined circular orbits is in- 
dependent of mass to first order, and is given by Eq. ()2.15|) . This result restricts the migratory 
behavior of EGPs, placing constraints on the initial semimajor axis range where a system can be 
dynamically excited. 2) For almost circular, arbitrary inclined orbits, we obtained analytical for- 
mulas for the secular evolution of some orbital parameters expressed in terms of the masses and 
initial semimajor axes only (Eqs. l4.12M.14"j 14.181 and I4.19j) . 3) Laplace-Lagrange secular theory 
fails to trace the eccentricity evolution two planets, both on initially circular orbits. 4) Libration 
widths computed with nonzero values of mass, eccentricity, and inclination for each planet (Eq. 
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I4.46|) appear to be viable tracers of the dynamical reach of resonances. 5) Encounter maps may be 
used to study the gross properties of the dynamics of EGP systems with a suitable scaling of the 
Hill parameter. 
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Tabic 1: Summary of the number of dynamically rigorous but stable systems affected by every 1st, 2nd, and 3rd 
order resonance in the semimajor axis range sampled in the simulations. The effect of each resonance is measured by 
its libration width, and twice that value, denoted respectively by lw and 2w. Initial relative inclinations are given 
by the leftmost column. 
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Initial Variables 


Final Variables 


Mapping Variables 


Other Variables 


■m ,m 1 ,m 2 






1 

_ mi+m 2 3 


a 1 ,a 2 


a r = ai - a2 


K=lal 


a 2 < a < ai 


ei,e 2 


r _ ei-e 2 
c r — e 


1 2 


e * = 




— «l-«2 
£ r — ^ 


J — 2 




Ai, A2 


A r = Al — A2 






"ij ^2 


= r^i — r^2 










h—e r cos w r 
k—e r sinov 





Table 2: Summary of variables used in the NLTP map. 
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Figure Captions 



Figure H Critical Hill Separation as a function of relative initial inclination when ay — 1. The 
dot-dashed, dashed, and dotted lines correspond to two planets with masses 1O~ 3 M , 1O~ 4 M , 
and 1O _5 M respectively. The solid lines indicate the locations of major first and second order 
resonances. 

Figure El Time to instability for initially unstable coplanar systems. Diamonds mark unstable 
systems, and bars mark stable systems, all of which were stopped after 2 x 10 6 . These time units 
correspond to years for an inner planet with a semimajor axis of 1 AU. 

Figure 01 Same as Fig. |21 for systems with initial relative inclinations of 0° (diamonds), 5° 
(triangles), 10° (squares), and 15° (crosses). Stable systems are not shown. 

Figure EJ Record of instability. Shown are four sets of four rows of bars, each set corresponding 
to a different initial relative inclination value, and each row corresponding to a different trial for 
given initial values of a 2 and I. Each bar represents an unstable simulation; absence of a bar 
indicates a stable simulation; an "x" represents simulations that the HNbody integrator failed to 
complete because it could not achieve the desired accuracy. The horizontal axis provides the initial 
planet separation, where the critical separation A crit is computed from Eq. (j2.15J) . The dot-dashed 
vertical lines represent global chaos limits (oc fj 2 ^ 7 ) conjectured by Wisdom (1980) for equal mass 
planets on circular planar orbits. The coefficient of /z 2 / 7 was fit for the 7 = 0° case, and applied to 
the other cases. The 5:3 resonance location is displayed with solid vertical lines. 

Figure 03 Continuation of Fig. 0]for initial relative inclination values of 20° to 35°. All first, 
second, and third order nominal resonances locations in the sampled regions are shown. The region 
of global chaos extends beyond (is less than) A crit /2 for most of the inclinations sampled here, and 
hence is not seen except around A cr i t /2 for / = 20°. 

Figure |H1 Same as Fig. 0] except that bars represent stable systems where significant migration 
(semimajor axis of either planet differing by at least 20% of its initial value) occurs. 

Figure [7| Continuation of Fig. with the locations of the relevant major resonances shown, 
including the fourth-order 7:3 resonance. 

Figure |H1 Final vs. initial semimajor axis of the outer planet for the stable systems with an 
initial relative inclination of 15° (upper panel), 25° (middle panel), and 30° (lower panel). 

Figure El Final values of eccentricity vs. semimajor axis for the inner planet (top) and outer 
planet (bottom) in the systems from the middle panel of Fig. |H1 

Figure E3 Comparison of outputted data from numerical simulations (crosses) and from Eqs. 
(j4.12|) - (j4.14|) (solid lines) for a high initial relative inclination (35°) system with an initial outer 
semimajor axis value only 0.0003 away from instability. 

Figure El Comparison of outputted data from numerical simulations (crosses) and from Eqs. 
(J4.18)) and 1)4.19)1 (solid lines) for a high inclination (35°) system with an initial outer semimajor 
axis value only 0.0003 away from instability. 
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Figure Variation of libration width over one dynamical inclination period for two different 
types of 5:3 resonances and 6 different relative initial inclinations. 

Figure Reproduction of Fig. 0] with results of the second-order nonplanar eccentric map 
superimposed. The additional, darker bars represent unstable systems from the map. 

Figure ITU Stability of planetary systems in the initial 02 interval [1, 2], for inclinations up to 
17.5°. Each dark bar represents an unstable system. 

Figure ITol Same as Fig. except for inclinations in the range [20°, 37.5°]. 
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